A C++ library with a set of Python scripts for medical image processing
2024-11-04
// Content of `math-tools/common_tools/convert_image/animaCreateImage.cxx`
#include <tclap/CmdLine.h>
#include <iostream>
#include <string>
#include <itkImage.h>
#include <itkVectorImage.h>
#include <itkImageRegionIterator.h>
#include <animaReadWriteFunctions.h>
int main(int argc, char **argv)
{
TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
TCLAP::ValueArg<std::string> geomArg("g","geometryfile","Geometry image",true,"","Geometry image",cmd);
TCLAP::ValueArg<std::string> outArg("o","outputfile","output image",true,"","output image",cmd);
TCLAP::ValueArg<unsigned int> vdimArg("v","vdim","Force vdim to this value",false,1,"Force vdim",cmd);
TCLAP::ValueArg<unsigned int> valueArg("b","buffervalue","Value to fill buffer",false,0,"Value to fill buffer",cmd);
TCLAP::SwitchArg vecArg("V","isvec","Input image is a vector / tensor image",cmd,false);
try
{
cmd.parse(argc,argv);
}
catch (TCLAP::ArgException& e)
{
std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
return EXIT_FAILURE;
}
using ImageType = itk::Image <double,3>;
using VectorImageType = itk::VectorImage <double,3>;
bool isVect = vecArg.isSet();
unsigned int fvdim = vdimArg.getValue();
if (fvdim > 1)
isVect = true;
if (!isVect)
{
ImageType::Pointer geomImage = anima::readImage <ImageType> (geomArg.getValue());
ImageType::Pointer resImage = ImageType::New();
resImage->Initialize();
resImage->SetRegions(geomImage->GetLargestPossibleRegion());
resImage->SetSpacing(geomImage->GetSpacing());
resImage->SetOrigin(geomImage->GetOrigin());
resImage->SetDirection(geomImage->GetDirection());
resImage->Allocate();
resImage->FillBuffer(valueArg.getValue());
anima::writeImage <ImageType> (outArg.getValue(),resImage);
}
else
{
VectorImageType::Pointer geomImage = anima::readImage <VectorImageType> (geomArg.getValue());
VectorImageType::Pointer resImage = VectorImageType::New();
resImage->Initialize();
resImage->SetRegions(geomImage->GetLargestPossibleRegion());
resImage->SetSpacing(geomImage->GetSpacing());
resImage->SetOrigin(geomImage->GetOrigin());
resImage->SetDirection(geomImage->GetDirection());
unsigned int vdim = geomImage->GetNumberOfComponentsPerPixel();
if (fvdim > 1)
vdim = fvdim;
resImage->SetNumberOfComponentsPerPixel(vdim);
resImage->Allocate();
itk::VariableLengthVector <double> tmpData;
tmpData.SetSize(vdim);
for (unsigned int i = 0;i < vdim;++i)
tmpData[i] = valueArg.getValue();
itk::ImageRegionIterator <VectorImageType> tmpIt(resImage,resImage->GetLargestPossibleRegion());
while(!tmpIt.IsAtEnd())
{
tmpIt.Set(tmpData);
++tmpIt;
}
anima::writeImage <VectorImageType> (outArg.getValue(),resImage);
}
return EXIT_SUCCESS;
}\[ \log \left( - \log(p) \right); \]
Tip
Get inspiration from (copy-paste) existing code and modify the minimum to suit your needs!
double_log_pvalue);CMakeLists.txt file by adding add_subdirectory(double_log_pvalue). This tells cmake to add the subdirectory for compilation;animaDoubleLogPValue.cxx) in the double_log_pvalue folder.Now compile your code.
Caution
If you try to compile this code, it will not work because your tool is missing its own CMakeLists.txt file to give instructions to cmake on how to compile our tool (i.e. which libraries should it link with, which headers, etc.).
Do I have to learn cmake as well?
Tip
As usual, let us adopt the minimum effort strategy and see if the existing CMakeLists.txt for animaCreateImage could fit our needs.
# Context of `math-tools/common_tools/create_image/CMakeLists.txt`
if(BUILD_TOOLS)
project(animaCreateImage)
## #############################################################################
## List Sources
## #############################################################################
list_source_files(${PROJECT_NAME}
${CMAKE_CURRENT_SOURCE_DIR}
)
## #############################################################################
## add executable
## #############################################################################
add_executable(${PROJECT_NAME}
${${PROJECT_NAME}_CFILES}
)
## #############################################################################
## Link
## #############################################################################
target_link_libraries(${PROJECT_NAME}
${ITKIO_LIBRARIES}
)
## #############################################################################
## install
## #############################################################################
set_exe_install_rules(${PROJECT_NAME})
endif()math-tools/common_tools/create_image/CMakeLists.txt in the double_log_pvalue folder;Note
Anima handles most of the compilation settings for you. You only need to care about the specific CMakeLists.txt for your tool but, even there, most of the time, copy-pasting an already existing one will be enough.
main() function;# Use a terminal to search for existing filters
# based on the one you want to use.
stamm-a@aymerics-mbp Anima % grep -r "animaMasked" *
diffusion/dti/dti_estimator/animaDTIEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
diffusion/dti/flip_tensors/animaFlipTensorImageFilter.h:#include <animaMaskedImageToImageFilter.h>
diffusion/dti/dti_non_central_chi_estimator/animaDTINonCentralChiEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
diffusion/mcm_estimator/animaMCMEstimatorImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/statistics/animaLocalPatchCovarianceDistanceImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/statistics/animaLocalPatchMeanDistanceImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/common/animaMaskedImageToImageFilter.hxx:#include "animaMaskedImageToImageFilter.h"
math-tools/common/animaMaskedImageToImageFilter.h:#include "animaMaskedImageToImageFilter.hxx"
math-tools/statistical_tests/animaCramersTestImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/statistical_tests/nlmeans_patient_to_group_comparison/animaNLMeansPatientToGroupComparisonImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/statistical_tests/patient_to_group_comparison/animaPatientToGroupComparisonImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/animaT1SERelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/gamma_mixture_t2_estimation/animaGammaMixtureT2RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/gmm_t2_estimation/animaGMMT2RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/animaT2EPGRelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/combined_relaxometry_estimation/animaCombinedRelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/multi_t2_estimation/animaMultiT2RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/animaT1RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/animaT2RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
registration/resamplers/animaOrientedModelBaseResampleImageFilter.h:#include <animaMaskedImageToImageFilter.h>
segmentation/staple/animaMultiThreadedSTAPLEImageFilter.h:#include <animaMaskedImageToImageFilter.h>
segmentation/tissues_em_classification/animaTissuesEMClassificationImageFilter.h:#include <animaMaskedImageToImageFilter.h>An image filter is a C++ class templated over the type of input images and the type of output images;
It is made of two files:
.h header file containing method declarations;.hxx file containing method implementations.Naming convention is animaToolNameImageFilter.h and animaToolNameImageFilter.hxx.
// Content of `math-tools/statistics/animaCramersTestImageFilter.h`
#pragma once
#include <iostream>
#include <animaMaskedImageToImageFilter.h>
#include <itkVectorImage.h>
#include <itkImage.h>
#include <vector>
namespace anima
{
template <class PixelScalarType>
class CramersTestImageFilter :
public anima::MaskedImageToImageFilter< itk::VectorImage <PixelScalarType, 3> , itk::Image <PixelScalarType, 3> >
{
public:
/** Standard class typedefs. */
using Self = CramersTestImageFilter<PixelScalarType>;
using TInputImage = itk::VectorImage <PixelScalarType, 3>;
using TOutputImage = itk::Image <PixelScalarType, 3>;
using Superclass = anima::MaskedImageToImageFilter< TInputImage, TOutputImage >;
using Pointer = itk::SmartPointer<Self>;
using ConstPointer = itk::SmartPointer<const Self>;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods) */
itkTypeMacro(CramersTestImageFilter, MaskedImageToImageFilter);
/** Image typedef support */
using InputImageType = TInputImage;
using OutputImageType = TOutputImage;
using InputPixelType = typename InputImageType::PixelType;
using InputImagePointer = typename InputImageType::Pointer;
using OutputImagePointer = typename OutputImageType::Pointer;
/** Superclass typedefs. */
using MaskImageType = typename Superclass::MaskImageType;
using MaskImagePointer = typename MaskImageType::Pointer;
using OutputImageRegionType = typename Superclass::OutputImageRegionType;
/** Set the number of samples to build the distribution */
itkSetMacro(NbSamples, unsigned long int);
itkGetMacro(NbSamples, unsigned long int);
void SetFirstGroupIndexes(std::vector<unsigned int> &group)
{
m_FirstGroup.clear();
m_FirstGroupSize = group.size();
for (unsigned int i = 0;i < m_FirstGroupSize;++i)
m_FirstGroup.push_back(group[i]);
}
void SetSecondGroupIndexes(std::vector<unsigned int> &group)
{
m_SecondGroup.clear();
m_SecondGroupSize = group.size();
for (unsigned int i = 0;i < m_SecondGroupSize;++i)
m_SecondGroup.push_back(group[i]);
}
void AddOutlierMask(MaskImageType *mask)
{
m_OutlierMasks.push_back(mask);
}
protected:
CramersTestImageFilter()
: Superclass()
{
m_NbSamples = 5000;
m_FirstGroup.clear();
m_SecondGroup.clear();
m_FirstGroupSize = 0;
m_SecondGroupSize = 0;
m_UseOutlierMasks = false;
m_OutlierMasks.clear();
}
virtual ~CramersTestImageFilter() {}
void BeforeThreadedGenerateData() ITK_OVERRIDE;
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE;
//! Bootstrap samples generator for group comparison
void GenerateBootStrapSamples();
//! Actually bootstraps a p-value from the computed distance matrix and group samples
double BootStrap(vnl_matrix <double> &groupDistMatrix, std::vector <double> &inlierWeights);
//! Computes the Cramers' statistic knowing the distance matrix and the two groups
double CramerStatistic(vnl_matrix <double> &grpDistMatrix, std::vector <double> &inlierWeights,
unsigned int index);
private:
ITK_DISALLOW_COPY_AND_ASSIGN(CramersTestImageFilter);
std::vector < unsigned int > m_FirstGroup, m_SecondGroup;
std::vector < std::vector < unsigned int > > m_SamplesFirstGroup, m_SamplesSecondGroup;
unsigned long int m_NbSamples;
unsigned int m_FirstGroupSize, m_SecondGroupSize;
std::vector <MaskImagePointer> m_OutlierMasks;
bool m_UseOutlierMasks;
};
} // end namespace anima
#include "animaCramersTestImageFilter.hxx"Modify public/protected functions & private attributes as needed to clarify your code.
// Content of `math-tools/statistics/animaCramersTestImageFilter.hxx`
#pragma once
#include "animaCramersTestImageFilter.h"
#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
#include <itkTimeProbe.h>
#include <itkProgressReporter.h>
#include <cmath>
#include <ctime>
#include <random>
namespace anima
{
template <class PixelScalarType>
void
CramersTestImageFilter<PixelScalarType>
::BeforeThreadedGenerateData ()
{
Superclass::BeforeThreadedGenerateData();
this->GetOutput()->FillBuffer(0);
// Checking consistency of the data and parameters
unsigned int nbInputs = this->GetNumberOfIndexedInputs();
if (nbInputs <= 1)
itkExceptionMacro("Error: Not enough inputs available... Exiting...");
if ((m_FirstGroupSize + m_SecondGroupSize) != nbInputs)
itkExceptionMacro("Groups data not clearly wrong... Exiting...");
m_SamplesFirstGroup.clear();
m_SamplesSecondGroup.clear();
this->GenerateBootStrapSamples();
m_UseOutlierMasks = (m_OutlierMasks.size() == nbInputs);
}
template <class PixelScalarType>
void
CramersTestImageFilter<PixelScalarType>
::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread)
{
using InIteratorType = itk::ImageRegionConstIterator < TInputImage >;
using OutRegionIteratorType = itk::ImageRegionIterator < OutputImageType >;
using MaskRegionIteratorType = itk::ImageRegionIterator < MaskImageType >;
OutRegionIteratorType outIterator(this->GetOutput(), outputRegionForThread);
MaskRegionIteratorType maskIterator (this->GetComputationMask(), outputRegionForThread);
std::vector <InIteratorType> inIterators(this->GetNumberOfIndexedInputs());
for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
inIterators[i] = InIteratorType(this->GetInput(i), outputRegionForThread);
std::vector <MaskRegionIteratorType> outlierMasksIterators(m_OutlierMasks.size());
if (m_UseOutlierMasks)
{
for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
outlierMasksIterators[i] = MaskRegionIteratorType(m_OutlierMasks[i],outputRegionForThread);
}
std::vector < InputPixelType > firstGroupData(m_FirstGroupSize), secondGroupData(m_SecondGroupSize);
vnl_matrix <double> cramerDistMatrix(this->GetNumberOfIndexedInputs(),this->GetNumberOfIndexedInputs(),0.0);
std::vector <double> inlierProbabilities(this->GetNumberOfIndexedInputs(),1);
unsigned int vectorSize = this->GetInput(0)->GetNumberOfComponentsPerPixel();
while (!outIterator.IsAtEnd())
{
if (maskIterator.Get() == 0)
{
outIterator.Set(0.0);
++outIterator;
++maskIterator;
for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
++inIterators[i];
if (m_UseOutlierMasks)
{
for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
++outlierMasksIterators[i];
}
continue;
}
// Voxel is in mask, so now let's go for it and compute the stats
for (unsigned int i = 0; i < m_FirstGroupSize;++i)
{
unsigned int indFirst = m_FirstGroup[i];
firstGroupData[i] = inIterators[indFirst].Get();
}
for (unsigned int i = 0; i < m_SecondGroupSize;++i)
{
unsigned int indSec = m_SecondGroup[i];
secondGroupData[i] = inIterators[indSec].Get();
}
// Now that data is grouped, compute the distance matrix
for (unsigned int i = 0; i < m_FirstGroupSize;++i)
{
for (unsigned int l = 0;l < m_SecondGroupSize;++l)
{
double dist = 0;
for (unsigned int j = 0;j < vectorSize;++j)
dist += (firstGroupData[i][j] - secondGroupData[l][j]) * (firstGroupData[i][j] - secondGroupData[l][j]);
cramerDistMatrix(i,m_FirstGroupSize + l) = sqrt(dist);
cramerDistMatrix(m_FirstGroupSize + l,i) = cramerDistMatrix(i,m_FirstGroupSize + l);
}
for (unsigned int l = i+1;l < m_FirstGroupSize;++l)
{
double dist = 0;
for (unsigned int j = 0;j < vectorSize;++j)
dist += (firstGroupData[i][j] - firstGroupData[l][j]) * (firstGroupData[i][j] - firstGroupData[l][j]);
cramerDistMatrix(i,l) = sqrt(dist);
cramerDistMatrix(l,i) = cramerDistMatrix(i,l);
}
}
for (unsigned int i = 0; i < m_SecondGroupSize;++i)
{
for (unsigned int l = i+1;l < m_SecondGroupSize;++l)
{
double dist = 0;
for (unsigned int j = 0;j < vectorSize;++j)
dist += (secondGroupData[i][j] - secondGroupData[l][j]) * (secondGroupData[i][j] - secondGroupData[l][j]);
cramerDistMatrix(m_FirstGroupSize + i,m_FirstGroupSize + l) = sqrt(dist);
cramerDistMatrix(m_FirstGroupSize + l,m_FirstGroupSize + i) = cramerDistMatrix(m_FirstGroupSize + i,m_FirstGroupSize + l);
}
}
if (m_UseOutlierMasks)
{
for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
inlierProbabilities[i] = (outlierMasksIterators[i].Get() == 0);
}
outIterator.Set(this->BootStrap(cramerDistMatrix,inlierProbabilities));
++outIterator;
++maskIterator;
for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
++inIterators[i];
if (m_UseOutlierMasks)
{
for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
++outlierMasksIterators[i];
}
this->IncrementNumberOfProcessedPoints();
}
}
template <class PixelScalarType>
void
CramersTestImageFilter<PixelScalarType>
::GenerateBootStrapSamples()
{
m_SamplesFirstGroup.resize(m_NbSamples + 1);
m_SamplesSecondGroup.resize(m_NbSamples + 1);
m_SamplesFirstGroup[0] = m_FirstGroup;
m_SamplesSecondGroup[0] = m_SecondGroup;
unsigned int nbFirstGroup = m_FirstGroup.size();
unsigned int nbSecondGroup = m_SecondGroup.size();
unsigned int nbSubjects = nbFirstGroup + nbSecondGroup;
unsigned int minNbGroup = std::min(nbFirstGroup,nbSecondGroup);
bool isFirstGroupMin = (nbFirstGroup <= nbSecondGroup);
std::vector <unsigned int> sampleGen, sampleGenOtherGroup;
std::mt19937 generator(time(0));
std::uniform_real_distribution<double> unifDistr(0.0, 1.0);
for (unsigned int i = 1;i <= m_NbSamples;++i)
{
sampleGen.clear();
sampleGenOtherGroup.clear();
unsigned int j = 0;
while (j < minNbGroup)
{
int tmpVal = std::min((int)(nbSubjects - 1),(int)floor(unifDistr(generator) * nbSubjects));
if (tmpVal < 0)
tmpVal = 0;
bool isAlreadyIndexed = false;
for (unsigned int k = 0; k < j;++k)
{
if (tmpVal == sampleGen[k])
{
isAlreadyIndexed = true;
break;
}
}
if (!isAlreadyIndexed)
{
sampleGen.push_back(tmpVal);
j++;
}
}
// New sample gen that is not already here... Add it to the result vectors
for (unsigned int j = 0;j < nbSubjects;++j)
{
bool isAlreadyIndexed = false;
for (unsigned int k = 0;k < minNbGroup;++k)
{
if (j == sampleGen[k])
{
isAlreadyIndexed = true;
break;
}
}
if (!isAlreadyIndexed)
sampleGenOtherGroup.push_back(j);
}
if (isFirstGroupMin)
{
m_SamplesFirstGroup[i] = sampleGen;
m_SamplesSecondGroup[i] = sampleGenOtherGroup;
}
else
{
m_SamplesFirstGroup[i] = sampleGenOtherGroup;
m_SamplesSecondGroup[i] = sampleGen;
}
}
}
template <class PixelScalarType>
double
CramersTestImageFilter<PixelScalarType>
::BootStrap(vnl_matrix <double> &groupDistMatrix, std::vector <double> &inlierWeights)
{
// First index is the true group separation (see GenerateBootStrap)
double dataVal = this->CramerStatistic(groupDistMatrix,inlierWeights,0);
std::vector <double> statsValues(m_NbSamples,0);
for (unsigned long int i = 0;i < m_NbSamples;++i)
statsValues[i] = this->CramerStatistic(groupDistMatrix,inlierWeights,i+1);
std::sort(statsValues.begin(),statsValues.end());
if (dataVal >= statsValues[statsValues.size() - 1])
return 0;
unsigned int position = 0;
while(position < statsValues.size())
{
if (dataVal <= statsValues[position])
break;
++position;
}
double resVal = 0;
if (position > 0)
{
resVal = m_NbSamples - position + (statsValues[position - 1] - dataVal) / (statsValues[position] - statsValues[position - 1]);
resVal /= m_NbSamples;
}
else
{
// Here statsValues[position - 1] doesn't exist, replacing with 0
resVal = m_NbSamples - position - dataVal / statsValues[position];
resVal /= m_NbSamples;
}
return resVal;
}
template <class PixelScalarType>
double
CramersTestImageFilter<PixelScalarType>
::CramerStatistic(vnl_matrix <double> &grpDistMatrix, std::vector <double> &inlierWeights,
unsigned int index)
{
unsigned int nbFirstGroup = m_SamplesFirstGroup[index].size();
unsigned int nbSecondGroup = m_SamplesSecondGroup[index].size();
double firstTerm = 0, secondTerm = 0, thirdTerm = 0;
double Z1 = 0, Z2 = 0, Z3 = 0;
for (unsigned int i = 0; i < nbFirstGroup;++i)
{
unsigned int indFirst = m_SamplesFirstGroup[index][i];
for (unsigned int j = 0;j < nbSecondGroup;++j)
{
double alpha = inlierWeights[indFirst]*inlierWeights[m_SamplesSecondGroup[index][j]];
firstTerm += alpha*grpDistMatrix(indFirst,m_SamplesSecondGroup[index][j]);
Z1 += alpha;
}
for (unsigned int j = i+1;j < nbFirstGroup;++j)
{
double alpha = 0.5*inlierWeights[indFirst]*inlierWeights[m_SamplesFirstGroup[index][j]];
secondTerm += 2*alpha*grpDistMatrix(indFirst,m_SamplesFirstGroup[index][j]);
Z2 += 2*alpha;
}
// Don't forget to add the diagonal term for Z2
Z2 += inlierWeights[indFirst]*inlierWeights[m_SamplesFirstGroup[index][i]];
}
firstTerm /= Z1;
secondTerm /= (2*Z2);
for (unsigned int i = 0; i < nbSecondGroup;++i)
{
int indSec = m_SamplesSecondGroup[index][i];
for (unsigned int j = i+1;j < nbSecondGroup;++j)
{
double alpha = inlierWeights[indSec]*inlierWeights[m_SamplesSecondGroup[index][j]];
thirdTerm += 2*alpha*grpDistMatrix(indSec,m_SamplesSecondGroup[index][j]);
Z3 += 2*alpha;
}
// Don't forget to add the diagonal term for Z3
Z3 += inlierWeights[indSec]*inlierWeights[m_SamplesSecondGroup[index][i]];
}
thirdTerm /= (2*Z3);
return firstTerm - secondTerm - thirdTerm;
}
} // end namespace animaWrite an image filter to handle the p-value transformation;
Modify animaDoubleLogPValue.cxx to call the filter instead of doing the computation directly in the main() function;
Add an optional argument to the parser for setting the number of cores to parallelize over.
Tip
Again, get inspiration from the existing code to understand how to modify the main() function. See for instance animaCramersTest.cxx.
Anima Hackathon #1 - Rennes - aymeric.stamm@cnrs.fr